Introduction

This document will contain all analyses and figures produced in R to be included in the manuscript. Note that these are not all the figures to be included, some were produced in PowerPoint and have been saved as image files.

# Raw data, endpoint data, and row objects for tcplfit2
load("../Data/Padilla_DNT60_lmr0_w_egid.Rdata")
load("~/Desktop/StephanieProjects/DNT60Analysis/pipelining/pipelined data/Padilla_DNT60_mc0.rda")
load("../Data/Padilla_DNT60_mc0_n.rda")
load("../Data/Padilla_DNT60_rows_n.rda")

# Results of curve fitting
load("../Data/Padilla_DNT60_tcplfits.rda")
load("../Data/Padilla_DNT60_tcpl_out.Rdata")
tcpl_out.dt <- as.data.table(DNT60_tcpl_out)

# Summarizing results of curve-fitting
load("../Data/Padilla_DNT60_Strong Effectors by Endpoint.Rdata")
load("../Data/Padilla_DNT60_BMC ranges.Rdata") # This data needs to be recompiled

# Results of literature review
load("../Data/Padilla_DNT60_Behavioral LOELs.Rdata")
load("../Data/Literature Review Summary.Rdata")
load("../Data/Literature Review Concentrations and Activity.Rdata")
load("../Data/Literature Review LOELs and Dev LOELs.Rdata") 

Example Locomotor Response Data

lmr0.egid[cpid=="Chlorpyrifos (ethyl)", cpid := "Chlorpyrifos"]

Example tcpl Curve-Fit

Results

Appearance of Endpoint Data

Figure 3: Raw Fluoxetine Vehicle Control Speed Data with Average Speed Endpoint Data Raw and Transformed

# Time Series for Fluoxetine vehicle control

# Units
unit.t = "2 min"
unit.mov = "cm"
unit.conc = paste0("\U03BC","M")
no.A = 10


# Identify movement columns of interest and extract
t.cols <- grep("vt", names(lmr0.egid), value = TRUE)
cols <- t.cols[(no.A+1):length(t.cols)]
A.cols <- t.cols[!(t.cols%in%cols)]
to.fit <- lmr0.egid[wllt=="v" & egid=="E6", -A.cols, with=FALSE]

# Melt data and designate Light and Dark time periods
to.fit_melt <- melt(to.fit[,-"egid"], id.vars = names(to.fit)[1:9], variable.name = "t", value.name = "speed")
to.fit_melt[, t := as.factor( gsub("vt","",t) )]
to.fit_melt[t%in%11:30, phase := "Light"]
to.fit_melt[t%in%31:50, phase := "Dark"]

# Create title, x- and y-axis titles, legend labels, and legend title
title.t <- paste0("Time (",unit.t,")")
label.y <- "Speed"
title.mean <- paste0(label.y, " (",unit.mov,"/",unit.t,")")

title.legend <- "Experimental Phase"

## Create x-axis breaks and labels
m <- 50
x.breaks <- c(11, 20, 30, 40, 50)
x.labels1 <- 2*seq(from=no.A,to=m,by=10)
x.labels2 <- x.labels1 - 2
x.labels <- paste(x.labels2, x.labels1, sep="-")

# Printing of graphic isn't working very well so insert new breaks to make space for edges
x.breaks <- c(5, x.breaks, 55)
x.labels <- c("a", x.labels, "b")

## plot
FlxTSeriesPlot <- ggplot(to.fit_melt) +
          geom_boxplot(aes(x=t,y=speed,fill=phase)) +
          scale_x_discrete(breaks = x.breaks, labels = x.labels) +
          scale_fill_manual(values=c("darkgrey","white")) +
          labs(x = title.t, y = title.mean, fill = "Experimental Phase") +
          theme_bw() +
          theme(text = element_text(size = 26))
# Raw endpoint data for Fluoxetine vehicle control

sample_endpoints <- lapply(mc0, function(table) {
  table.egid <- gabi::data_egids(table)
  table.egid[egid=="E6"]
})

sample_endpoints[c("AUC_L","AUC_D","AUC_T")] = NULL
for (i in 1:13) {
  sample_endpoints[[i]][, endp := names(sample_endpoints)[i]]
}

endp.data <- do.call('rbind', sample_endpoints)
max.conc <- max( endp.data[cpid=="Fluoxetine", conc] )

raw_stats <- endp.data[wllt=="v" | (conc==max.conc & cpid=="Fluoxetine")][endp == "avgS_L"][, .(mean=mean(rval,na.rm=T),median=median(rval,na.rm=T)), by=.(conc)]
raw_stats_l <- data.table::melt(raw_stats, id.vars = "conc", variable.name = "stat")

FlxAvgS_LPlot <- ggplot(endp.data[wllt=="v" | (conc==max.conc & cpid=="Fluoxetine")][endp == "avgS_L"]) +
  geom_density(aes(x=rval, fill=as.factor(conc), color = as.factor(conc)), alpha = 0.5) +
  theme_bw() +
  theme(text = element_text(size = 26)) +
  scale_color_manual(values = colors[c(1,5)]) +
  scale_fill_manual(values = colors[c(1,5)]) +
  labs(color = paste0("Concentration (", unit.conc, ")"),
       fill = paste0("Concentration (", unit.conc, ")"),
       linetype = "",
       x = "Raw Endpoint Value",
       y = "Density")

FlxAvgS_LPlot <- FlxAvgS_LPlot + 
                  geom_vline(data = raw_stats_l, aes(xintercept=value, linetype=stat), color = colors[c(1,5,1,5)], size=1) +
                  scale_linetype_discrete(labels = c("Mean","Median")) +
                  scale_y_continuous(labels = c("0.0","0.1","0.2","0.3","")) +
                  labs(color = paste0("Concentration (", unit.conc, ")"),
                       fill = paste0("Concentration (", unit.conc, ")"),
                       linetype = "",
                       x = "Raw Endpoint Value",
                       y = NULL)

legend <- cowplot::get_legend(FlxAvgS_LPlot)

FlxAvgS_LPlot <- FlxAvgS_LPlot + theme(legend.position="none")
sample_endpoints_n <- lapply(mc0_n, function(list) {
  table <- list[[1]]
  table.egid <- gabi::data_egids(table)
  table.egid[egid=="E6"]
})

sample_endpoints_n[c("AUC_L","AUC_D","AUC_T")] = NULL
for (i in 1:13) {
  sample_endpoints_n[[i]][, endp := names(sample_endpoints_n)[i]]
}

endp.data_n <- do.call('rbind', sample_endpoints_n)

trans_stats <- endp.data_n[wllt=="v" | (conc==max.conc & cpid=="Fluoxetine")][endp == "avgS_L"][, .(mean=mean(rval,na.rm=T),median=median(rval,na.rm=T)), by=.(conc)]
trans_stats_l <- data.table::melt(trans_stats, id.vars = "conc", variable.name = "stat")

FlxAvgS_LPlot_n <-ggplot(endp.data_n[wllt=="v" | (conc==max.conc & cpid=="Fluoxetine")][endp == "avgS_L"]) +
  geom_density(aes(x=rval, fill=as.factor(conc), color = as.factor(conc)), alpha = 0.5) +
  theme_bw() +
  theme(text = element_text(size = 26)) +
  scale_color_manual(values = colors[c(1,5)]) +
  scale_fill_manual(values = colors[c(1,5)]) +
  labs(color = paste0("Concentration (", unit.conc, ")"),
       fill = paste0("Concentration (", unit.conc, ")"),
       linetype = "",
       x = "Transformed Endpoint Value",
       y = NULL)

FlxAvgS_LPlot_n <- FlxAvgS_LPlot_n + 
                    geom_vline(data = trans_stats_l, aes(xintercept=value, linetype=stat), color = colors[c(1,5,1,5)], size=1) +
                    scale_linetype_discrete(labels = c("Mean","Median")) +
                    labs(color = paste0("Concentration (", unit.conc, ")"),
                         fill = paste0("Concentration (", unit.conc, ")"),
                         linetype = "",
                         x = "Transformed Endpoint Value",
                         y = NULL)

FlxAvgS_LPlot_n <- FlxAvgS_LPlot_n + theme(legend.position="none")
FlxEndpPlots <- cowplot::plot_grid(FlxAvgS_LPlot, FlxAvgS_LPlot_n,
                   labels = list("B","C"), label_size = 24)
FlxEndpPlots1 <- cowplot::plot_grid(FlxEndpPlots, legend, rel_widths = c(3, .4))
FlxPlots <- cowplot:::plot_grid(FlxTSeriesPlot, FlxEndpPlots1,
                    labels = list("A"), label_size = 24,
                    ncol = 1)
FlxPlots

tiff(filename = "test.tiff", height = 5, width = 5, units = "in", res=300)
FlxPlots
dev.off()
## png 
##   2

Supplemental Figure 1: Fluoxetine Vehicle Control and High Concentration Group Raw Endpoint Data

# Raw endpoint data for Fluoxetine vehicle control
raw_stats1 <- endp.data[wllt=="v" | (conc==max.conc & cpid=="Fluoxetine")][, .(mean=mean(rval,na.rm=T),median=median(rval,na.rm=T)), by=.(endp,conc)]
raw_stats1_l <- data.table::melt(raw_stats1, id.vars = c("endp","conc"), variable.name = "stat")

raw_plot_facet <- ggplot(endp.data[wllt=="v" | (conc==max.conc & cpid=="Fluoxetine")]) +
                    geom_density(aes(x=rval, fill=as.factor(conc), color = as.factor(conc)), alpha = 0.5) +
                    theme_bw() +
                    theme(text = element_text(size = 16)) +
                    facet_wrap(~endp, scales = "free") +
                    scale_color_manual(values = colors[c(1,5)]) +
                    scale_fill_manual(values = colors[c(1,5)])

raw_plot_facet <- raw_plot_facet + 
                    geom_vline(data = raw_stats1_l, aes(xintercept=value, linetype=stat), color = rep(colors[c(1,5)],26), size=1) +
                                scale_linetype_discrete(labels = c("Mean","Median")) +
                                labs(title = "Raw Endpoint Data for Fluoxetine Vehicle Control and High Concentration Group",
                                       color = paste0("Concentration (", unit.conc, ")"),
                                       fill = paste0("Concentration (", unit.conc, ")"),
                                       x = "Raw Endpoint Value",
                                       y = "Density")
raw_plot_facet

Supplemental Table 1: Calculate means and variance for endpoints and stick into table.

sample.endp.stats <- endp.data[wllt=="v", .(mean=mean(rval,na.rm=T), variance=stats::var(rval,na.rm=T)), by=.(endp)]
sample.endp.stats[, `:=` (mean = signif(mean,digits=3), variance = signif(variance,digits=3))]
DT::datatable(sample.endp.stats, colnames=c("Endpoint Abbreviation","Mean","Variance"),
              caption = "Raw Endpoint Data Statistics for Fluoxetine Vehicle Control")

Supplemental Table 2: Box-Cox power transformation parameters.

bxcx.params <- lapply(mc0_n, function(data) {
  as.data.table( data[c(2,3)] )
})

bxcx.params.dt <- do.call('rbind', bxcx.params)
bxcx.params.dt[, endp := names(bxcx.params)]

datatable(bxcx.params.dt[, .(endp,lam.hat,shift)], colnames=c("Endpoint Abbreviation","Lambda","Shift"),
              caption = "Box-Cox Power Transformation Parameters for Endpoints")

Supplemental Figure 2: Fluoxetine Vehicle Control and High Concentration Group Transformed Endpoint Data

trans_stats1 <- endp.data_n[wllt=="v" | (conc==max.conc & cpid=="Fluoxetine")][, .(mean=mean(rval,na.rm=T),median=median(rval,na.rm=T)), by=.(endp,conc)]
trans_stats1_l <- data.table::melt(trans_stats1, id.vars = c("endp","conc"), variable.name = "stat")

trans_plot_facet <- ggplot(endp.data_n[wllt=="v" | (conc==max.conc & cpid=="Fluoxetine")]) +
                      geom_density(aes(x=rval, fill=as.factor(conc), color = as.factor(conc)), alpha = 0.5) +
                      theme_bw() +
                      theme(text = element_text(size = 16)) +
                      facet_wrap(~endp, scales = "free") +
                      scale_color_manual(values = colors[c(1,5)]) +
                      scale_fill_manual(values = colors[c(1,5)])

trans_plot_facet <- trans_plot_facet + 
                      geom_vline(data = trans_stats1_l, aes(xintercept=value, linetype=stat), color = rep(colors[c(1,5)],26), size=1) +
                                scale_linetype_discrete(labels = c("Mean","Median")) +
                                labs(title = "Raw Endpoint Data for Fluoxetine Vehicle Control and High Concentration Group",
                                       color = paste0("Concentration (", unit.conc, ")"),
                                       fill = paste0("Concentration (", unit.conc, ")"),
                                       x = "Raw Endpoint Value",
                                       y = "Density")

trans_plot_facet

Applying tcplFit2

How many chemical-endpoint pair were active? How many chemicals had at least one active endpoint?

# Endpoints to exclude
exclude <- c("AUC_L","AUC_D","AUC_T")

# Number of active chemical-endpoint pairs
tcpl_out.dt[hitcall>0.8 & !endp%in%exclude, .N]
## [1] 57
# Number of active chemicals
length( tcpl_out.dt[hitcall>0.8 & !endp%in%exclude, unique(name)] )
## [1] 22

Evaluating popular fits.

# Popular fits for all endpoint data
tcpl_out.dt[!endp%in%exclude, .N, by = .(fit_method)]
##    fit_method   N
## 1:       exp4 272
## 2:      poly1 412
## 3:       exp5  16
## 4:       hill   5
## 5:       gnls  47
## 6:        pow  38
## 7:       exp2   3
# Popular fits for active curve-fits
tcpl_out.dt[!endp%in%exclude & hitcall>0.8, .N, by = .(fit_method)]
##    fit_method  N
## 1:       exp4 11
## 2:      poly1 22
## 3:       exp5  5
## 4:       gnls  6
## 5:        pow 10
## 6:       exp2  1
## 7:       hill  2

Evaluate the the occurrence of developmental toxicity for chemicals fit optimally by the poly1 function.

# Isolate chemicals with active endpoints optimally fitted by poly1
chemicals <- tcpl_out.dt[hitcall>0.8 & fit_method=="poly1" & !endp%in%c("AUC_L","AUC_D","AUC_T"), unique(name)]

# Is developmental toxicity associated with these chemicals?
litReviewConc[, dev_LOEL := min(.SD[devtox==1,conc]), by=.(cpid)]
litReviewConc[is.infinite(dev_LOEL), dev_LOEL := NA]
litReviewConc[cpid%in%chemicals, .(dev_LOEL = unique(dev_LOEL)), by = .(cpid)]
##                                        cpid dev_LOEL
##  1:                          5-Fluorouracil       NA
##  2:                             Amphetamine    120.0
##  3:                       Bisphenol A (BPA)    120.0
##  4:                            Chlorpyrifos     12.0
##  5:                                Dieldrin      0.4
##  6:                     Diethylstilbesterol     10.0
##  7:                              D-sorbitol       NA
##  8:                              Fluoxetine     12.0
##  9:                              Heptachlor      4.0
## 10:                              Loperamide    120.0
## 11: Polybrominated diphenyl ether (PBDE)-47     13.0
## 12:                            Tebuconazole     40.0

Activity in Terms of Endpoints

Table 2: Number of chemicals active in each endpoint.

count.hits.endp <- tcpl_out.dt[hitcall>0.8 & !endp%in%exclude, unique(name), by = .(endp)][, .N, by = .(endp)]

datatable(count.hits.endp[order(-N)], colnames=c("Endpoint Abbreviation","Number of Active Chemicals"), caption = "Number of Chemicals Active in Each Endpoint",
          rownames = FALSE)

Figure 4: Number of chemicals active for each endpoint category and experimental phase.

tcpl_out.actives <- tcpl_out.dt[hitcall>0.8 & !endp%in%c("AUC_L","AUC_D","AUC_T","AUC_r")]

tcpl_out.actives[endp%in%c("avgS_L","avgS_D","avgS_T"), cat1 := "Average Speed"]
tcpl_out.actives[endp%in%c("hbt1_L","hbt1_D","hbt2_L","hbt2_D"), cat1 := "Habituation"]
tcpl_out.actives[endp%in%c("RoA_L","RoA_D"), cat1 := "Range of Activity"]
tcpl_out.actives[endp%in%c("strtlA","strtlAavg","strtlF"), cat1 := "Startle Response"]

tcpl_out.actives[endp%in%c("avgS_L","hbt1_L","hbt2_L","RoA_L"), cat2 := "Light"]
tcpl_out.actives[endp%in%c("avgS_D","hbt1_D","hbt2_D","RoA_D"), cat2 := "Dark"]
tcpl_out.actives[endp%in%c("strtlA","strtlAavg","strtlF"), cat2 := "Transition"]
tcpl_out.actives[endp%in%c("avgS_T"), cat2 := "Total"]

actv.chem.cat1 <- tcpl_out.actives[, .(cat1 = factor(unique(cat1),levels=c("Average Speed","Habituation",
                                                                           "Range of Activity","Startle Response"))),
                                   by=.(name)]
actv.chem.cat2 <- tcpl_out.actives[, .(cat2 = factor(unique(cat2),levels=c("Light","Transition","Dark","Total"))),
                                   by=.(name)]

activesCatCntPlot <- ggplot(actv.chem.cat1, aes(x=cat1, color=cat1, fill=cat1)) +
  geom_bar() +
  stat_count(geom = "text", colour = "black", size = 3.5,
             aes(label = paste(..count..,"22",sep="/")), position=position_stack(vjust=1.05)) +
  theme_bw() +
  theme(text = element_text(size=14)) + 
  labs(x="Endpoint Categories",
       y="Number of Chemicals Active in Category",
       color="Endpoint Category",
       fill="Endpoint Category") +
  scale_y_continuous(breaks=c(2,4,6,8,10)) +
  scale_x_discrete(labels=c("Average Speed","Habituation","Range of Activity","Startle Response")) +
  scale_color_manual(values=viridis::viridis(4),
                     labels=c("Average Speed","Habituation","Range of Activity","Startle Response")) +
  scale_fill_manual(values=viridis::viridis(4),
                    labels=c("Average Speed","Habituation","Range of Activity","Startle Response"))


activesPhsCntPlot <- ggplot(actv.chem.cat2, aes(x=cat2, color=cat2, fill=cat2)) +
                            geom_bar() +
  stat_count(geom = "text", colour = "black", size = 3.5,
             aes(label = paste(..count..,"22",sep="/")), position=position_stack(vjust=1.05)) +
                            theme_bw() +
                            theme(text = element_text(size=14)) + 
                            labs(x="Experimental Phase",
                                 y="Number of Chemicals Active in Phase",
                                 color="Experimental Phase",
                                 fill="Experimental Phase") +
                                 scale_y_continuous(breaks=c(2,4,6,8,10,12)) +
                                 scale_color_manual(values=viridis::viridis(4)) +
                                 scale_fill_manual(values=viridis::viridis(4))

cowplot:::plot_grid(activesCatCntPlot, activesPhsCntPlot,
                    labels = "AUTO", label_size = 24,
                    ncol = 2)

Supplemtal Figure 3: Combinations of Endpoint Category and Experimental Phase Activity

endp.lists <- split(actv.chem.cat1[,cat1], actv.chem.cat1[,name])
to.plot <- data.table(name = names(endp.lists), list = endp.lists)

phase.lists <- split(actv.chem.cat2[,cat2], actv.chem.cat2[,name])
to.plot1 <- data.table(name = names(phase.lists), list = phase.lists)

plot <- ggplot(to.plot, aes(x=list)) +
          geom_bar() +
          theme_bw() +
          theme(text = element_text(size=14)) +
          scale_x_upset(reverse = TRUE) +
          labs(x = "Endpoint Category Combinations",
               y = "Number of Chemicals with Combination",
               title = "Combinations of Endpoint Activity by Endpoint Category")

plot1 <- ggplot(to.plot1, aes(x=list)) +
          geom_bar() +
          theme_bw() +
          theme(text = element_text(size=14)) +
          ggupset::scale_x_upset(reverse = TRUE) +
          labs(x = "Experimental Phase Combinations",
               y = "Number of Chemicals with Combination",
               title = "Combinations of Endpoint Activity by Experimental Phase")

cowplot::plot_grid(plot, plot1, rel_widths = c(9,8),
                   labels = "AUTO", label_size = 24)

Directionality of endpoint activity was taken from BMC heatmap which will be presented in BMC section

Results in Terms of Chemicals

Chemicals that were active in endpoints other than Average Speed endpoints.

active.avgS <- tcpl_out.dt[hitcall>0.8 & endp%in%c("AUC_L","AUC_D","AUC_T"), unique(name)]
tcpl_out.dt[hitcall>0.8 & !name%in%active.avgS & !endp%in%exclude, unique(name)]
##  [1] "5,5-Diphenylhydantoin"  "5-Fluorouracil"         "Acrylamide"            
##  [4] "Bis(tributyltin) oxide" "Chlorpyrifos"           "Diethylstilbesterol"   
##  [7] "Dieldrin"               "D-sorbitol"             "Heptachlor epoxide"    
## [10] "Phenobarbital"          "Tebuconazole"           "Triethyltin"

Linear Regression of Chlorpyrifos Average Speed in Dark

cpf_row <- rows_n[["Chlorpyrifos"]][["avgS_D"]]

cpf_dt <- data.frame(conc = c(cpf_row$conc,rep(0,length(cpf_row$bresp))),
                     resp = c(cpf_row$resp,cpf_row$bresp))
cpf_dt$conc = cpf_dt$conc - mean(cpf_dt$conc)

# Regress responses with a linear model
cpf_lm <- lm(resp ~ conc, data = cpf_dt)
confint(cpf_lm)
##                   2.5 %     97.5 %
## (Intercept) -0.02297693 0.04307083
## conc        -0.04239414 0.01099271

Figure 5: Behavior profiles of Chlorpyrifos with active endpoints.

chemical <- "Chlorpyrifos"

# Units
unit.t = "min"
unit.mov = "cm"
unit.conc = paste0("\U03BC","M")
prsp = "SA"
no.A = 10

# Extract chemical data
group <- unique(lmr0.egid[cpid == chemical, egid])

## Identify movement columns of interest
t.cols <- grep("vt", names(lmr0.egid), value = TRUE)
cols <- t.cols[(no.A+1):length(t.cols)]
A.cols <- t.cols[!(t.cols%in%cols)]

## extract data to be plotted, exclude acclimation
to.fit <- lmr0.egid[cpid==chemical | (wllt=="v" & egid==group), -A.cols, with=FALSE]

# create appropriate axes titles for plots
label.y <- "Speed"

# Format data for plotting

## calculate mean and 50% CIs for each vector column by concentration group, excluding concentration
exclude.A <- t.cols[!(t.cols%in%A.cols)]
means <- to.fit[, lapply(.SD, function(col) mean(col,na.rm=T)),
                .SDcols = exclude.A,
                by = conc]

## calculate CI's for transformed values then transform back
shift <- 1
logCIs <- to.fit[, lapply(.SD, function(x) log10(x+shift)), .SDcols=exclude.A, by=conc][
  , lapply(.SD, function(x) t.test(x,conf.level=0.50)$conf.int), .SDcols=exclude.A, by=conc]
CIs <- logCIs[, lapply(.SD, function(x) (10^x)-shift), by=conc][,lapply(.SD, function(col) abs(diff(col))/2), .SDcols=exclude.A, by=conc]

## elongate means and CIs data, and join
means_long <- data.table::melt(means, id.vars = "conc", variable.name = "t", value.name = "mean")
means_long[, t := sub("vt","",t)]
CIs_long <- data.table::melt(CIs, id.vars = "conc", variable.name = "t", value.name = "CI")
CIs_long[, t := sub("vt","",t)]
stats <- means_long[CIs_long, on = c("conc","t")][, conc := as.factor(conc)]
stats[, t := as.numeric(t)]

# create standard error of mean estimates by time period and plot as ribbons or error bars

# plot time-series data

## create title, x- and y-axis titles, legend labels, and legend title
title <- paste0("Sample Averaged Time-Series for ", chemical)
title.t <- paste0("Time (",unit.t,")")
title.mean <- paste0("Mean ", label.y, " (",unit.mov,"/",unit.t,")")
conc.n <- to.fit[wllq==1, .N, by=.(conc)][order(conc)]
legend.labels <- paste0(conc.n$conc, ", n=", conc.n$N)

title.legend <- paste0("Concentration (", unit.conc, ")")

## get better colors for plotting
N <- length(unique(to.fit[,conc]))
colors <- viridis::viridis(N)

## create x-axis breaks and labels
m <- as.integer(max(means_long[,t]))
x.breaks <- seq(from=no.A,to=m,by=10)
x.labels1 <- 2*seq(from=no.A,to=m,by=10)
x.labels2 <- x.labels1 - 2
x.labels <- paste(x.labels2, x.labels1, sep="-")

## plot
SAplot_CPF <- ggplot() +
          geom_point(data = stats, aes(x=t, y=mean, color=as.factor(conc))) +
          geom_line(data = stats, aes(x=t, y=mean, color=conc, group=conc)) +
          scale_x_continuous(breaks = x.breaks, labels = x.labels) +
          scale_color_manual(values = colors, labels=legend.labels) +
          geom_ribbon(data = stats,
                      aes(x=t, ymax=mean+CI, ymin=mean-CI, group=conc, fill=conc),
                      alpha = 0.25) +
          geom_rect(aes(xmin=10,xmax=30,ymin=-1.25,ymax=-0.5) ,fill="white", color="black") +
          geom_rect(aes(xmin=30,xmax=50,ymin=-1.25,ymax=-0.5) ,fill="black", color="black") +
          annotate("text", x=c(20,40), y=rep(-0.875,2), color=c("black","white"), label=c("Light Phase","Dark Phase"), size=5) +
          scale_fill_manual(values = colors, labels=legend.labels) +
          labs(title = title, subtitle = "Acclimation Period Excluded: 50% Confidence Bands",
               x = title.t, y = title.mean, color = title.legend) +
          guides(fill = "none") +
          theme_bw() +
          theme(text = element_text(size = 18))
row1_CPF <- rows_n[[chemical]][["strtlAavg"]]
row2_CPF <- rows_n[[chemical]][["hbt1_D"]]

fit1_CPF <- concRespCoreZR(row1_CPF, do.plot = TRUE, verbose.plot = FALSE)[["plot"]]
fit2_CPF <- concRespCoreZR(row2_CPF, do.plot = TRUE, verbose.plot = FALSE)[["plot"]]

# fit1_CPF$theme$text$size = 4
# fit2_CPF$theme$text$size = 4
# 
# fit1_CPF$heme$line$size = 0.15
# grobs <- list(SAplot_CPF, fit1_CPF, fit2_CPF)
# layout <- matrix(c(1,1,1,1,2,3), nrow=3, ncol=2, byrow=TRUE)
# 
# gridExtra::grid.arrange(grobs = grobs, layout_matrix = layout)
p1 <- cowplot::plot_grid(fit1_CPF, fit2_CPF, nrow=1)
cowplot::plot_grid(SAplot_CPF, p1, nrow=2)

BMC Heatmap to be displayed later.

ac <- c("avgS_L", "hbt1_L", "hbt2_L","RoA_L",
        "strtlA","strtlAavg", "strtlF",
        "avgS_D","hbt1_D", "hbt2_D", "RoA_D",
        "avgS_T","AUC_r")

## Gather BMC's
bmd <- lapply(tcplfits_n, function(chm) unlist(lapply(chm, function(fit) {
  hitcall <- fit[["summary"]]$hitcall
  bmd <- fit[["summary"]]$bmd
  if (!(hitcall>0.8 & !is.na(bmd) & bmd != 0)) {
    bmd <- 10000
  }
  bmd
})))
bmd <- do.call("rbind", bmd)
bmd <- log10(bmd[,colnames(bmd)%in%ac])
# bmd <- bmd[,-c(1,9,14)]

# Identify active chemicals
actives <- rownames(bmd)[apply(bmd, 1, function(row) any(row < 4))]

# Isolate data of interest
to.fit <- bmd

# PBDE-47 name is too long
rownames(bmd)[50] <- "PBDE-47"
actives[19] <- "PBDE-47"
# Create a matrix specifying if an up, down, up and down, or down and up arrow should be printed in cells.
to.fit <- lapply(tcplfits_n, function(chm) chm[-c(1,9,14)])
layer.mat <- lapply(to.fit, function(chm) unlist(lapply(chm, function(fit) {
                hitcall <- fit[["summary"]]$hitcall
                fit_method <- fit[["summary"]]$fit_method
                dir <- sign(fit[["summary"]]$top)
                if (hitcall>0.8 & fit_method=="gnls" & dir==1) {
                  layer <- 0
                } else if (hitcall>0.8 & fit_method=="gnls" & dir==-1) {
                  layer <- 1
                } else if (hitcall>0.8 & dir==1) {
                  layer <- 2
                } else if (hitcall>0.8 & dir==-1) {
                  layer <- 3
                } else layer <- 4
              })))
layer.mat <- do.call("rbind", layer.mat)
rownames(layer.mat)[50] <- "PBDE-47" # Change PBDE-47 name to match BMD matrix
layer.mat1 <- layer.mat[row.names(layer.mat) %in% actives,]
# Column labels
col_labels <- c(expression("Average Speed in Light"^1), "Habiutation 1 in Light", "Habituation 2 in Light", expression("Range of Activity in Light"),
                "Startle Acceleration", "Startle Relative to Avg. Speed in Light", "Startle Fold-Change",
                expression("Average Speed in Dark"^1), "Habituation 1 in Dark", "Habituation 2 in Dark", expression("Range of Activity in Dark"),
                "Average Speed in Both Phases", expression("AUC in Dark / AUC in Light Ratio"^2)) # Superscripts notate references in poster
# Legends

# Custom heat legend.
f2 = circlize::colorRamp2(seq(min(bmd), max(bmd), length = 8), rev(viridis(8)), space = "sRGB")
heat_lgd = Legend(col_fun = f2,
                  title = paste0("BMC log(","\U03BC","M)"),
                  title_position = "lefttop",
                  legend_width = unit(4,"cm"),
                  direction = "horizontal")
# Column annotation by phase legend.
ann_lgd = Legend(labels = c("Light","Transition","Dark","Light+Dark"),
                 title = "Phase",
                 title_position = "leftcenter",
                 labels_gp = gpar(fontsize=8),
                 title_gp = gpar(fontsize=8),
                 legend_gp = grid::gpar(fill = c("white","grey","black","red")),
                 border = TRUE,
                 nrow = 1,
                 column_gap = unit(5, 'mm'))
# Legend for cell signal arrows.
dir_lgd = Legend(labels = c(paste("\U2191","Gain"),
                            paste("\U2193","Loss"),
                            paste("\U21C5","GainLoss"),
                            paste("\U21F5","LossGain")),
                 title = "Signal Direction",
                 labels_gp = gpar(fontsize=8),
                 title_gp = gpar(fontsize=8),
                 title_position = "leftcenter",
                 nrow = 1,
                 column_gap = unit(0, 'mm')) # Will produce warnings, don't worry.
lgd_list <- packLegend(ann_lgd, dir_lgd)

# Create column annotation indicating the phase of the LMR that is described.
column_ha <- ComplexHeatmap::HeatmapAnnotation(Phase = factor(c(rep("Light",4), rep("Transition",3),rep("Dark",4),rep("Light+Dark",2)),
                                                        levels=c("Light","Transition","Dark","Light+Dark")),
                                               border = TRUE,
                                               simple_anno_size = unit(0.25, 'cm'),
                                               col = list(Phase=c("Light"="white",
                                                                  "Transition"="grey",
                                                                  "Dark"="black",
                                                                  "Light+Dark"="red")),
                                               annotation_legend_param = list(nrow = 1),
                                               show_annotation_name = FALSE,
                                               show_legend = FALSE)

# Create function to add arrows indicating signal direction in cells.
cell_fun <- function(j, i, x, y, width, height, fill) {
  if (layer.mat1[i,j] == 0) {
    grid.text("\U21C5", x, y, gp=gpar(fontsize=6))
  } else if (layer.mat1[i,j] == 1) {
    grid.text("\U21F5", x, y, gp=gpar(fontsize=6))
  } else if (layer.mat1[i,j] == 2) {
    grid.text("\U2191", x, y, gp=gpar(fontsize=6))
  } else if (layer.mat1[i,j] == 3) {
    grid.text("\U2193", x, y, gp=gpar(fontsize=6))
  }
}
# Create main heat map.
htlist <- ComplexHeatmap::Heatmap(bmd[row.names(bmd) %in% actives,],

                        # Specify some parameters for heat legend.
                        name = paste0("BMC log(","\U03BC","M)"),
                        col = f2,
                        border_gp = grid::gpar(col="black",lwd=1),
                        rect_gp=grid::gpar(col="grey"),
                        show_heatmap_legend = TRUE,
                        heatmap_legend_param = list(legend_height = unit(3.5,"cm"),
                                                    direction = "vertical",
                                                    title_gp = grid::gpar(fontsize=10),
                                                    labels_gp = grid::gpar(fontsize=10)),

                        # Heatmap width and height
                        # width = unit(16,"in"),
                        # height = unit(18,"in"),

                        # Column label parameters
                        column_labels = col_labels, column_names_rot = 45,

                        # Split columns by phase.
                        column_split = factor(c(rep("Light",4), rep("Transition",3),rep("Dark",4),rep("Light+Dark",2)),
                                              levels=c("Light","Transition","Dark","Light+Dark")),

                        # Specify some column parameters.
                        column_title = NULL,
                        top_annotation = column_ha,

                        # Add signal direction arrows.
                        cell_fun = cell_fun,

                        # Specify some parameters row dendrogram aesthetics and row labels.
                        row_title_side = "right",
                        row_title_rot = 0,
                        row_split = 5,
                        row_dend_side = "right",
                        row_names_side = "left",

                        # Clustering parameters.
                        cluster_columns = FALSE,
                        cluster_rows = TRUE,
                        clustering_distance_rows = "pearson",

                        # Font sizes
                        row_names_gp = grid::gpar(fontsize=8),
                        column_names_gp = grid::gpar(fontsize=10),
                        row_title_gp = grid::gpar(fontsize=10),
                        column_title_gp = grid::gpar(fontsize=10)
                        )

Supplemental Figure 4: Activity in endpoints other than Average Speed endpoints.

isolate <- actives[order(actives)]
isolate <- c(isolate[c(1,2,6,8,21)], isolate[c(10,12,13,20,22)])

layer.mat2 <- layer.mat[isolate,]

cell_fun1 <- function(j, i, x, y, width, height, fill) {
  if (layer.mat2[i,j] == 0) {
    grid.text("\U21C5", x, y, gp=gpar(fontsize=12))
  } else if (layer.mat2[i,j] == 1) {
    grid.text("\U21F5", x, y, gp=gpar(fontsize=12))
  } else if (layer.mat2[i,j] == 2) {
    grid.text("\U2191", x, y, gp=gpar(fontsize=12))
  } else if (layer.mat2[i,j] == 3) {
    grid.text("\U2193", x, y, gp=gpar(fontsize=12))
  }
}
htlist2 <- ComplexHeatmap::Heatmap(bmd[isolate,],

                        # Specify some parameters for heat legend.
                        name = paste0("BMC log(","\U03BC","M)"),
                        col = f2,
                        border_gp = grid::gpar(col="black",lwd=1),
                        rect_gp=grid::gpar(col="grey"),
                        show_heatmap_legend = TRUE,
                        heatmap_legend_param = list(legend_height = unit(3.5,"cm"),
                                                    direction = "vertical",
                                                    title_gp = grid::gpar(fontsize=10),
                                                    labels_gp = grid::gpar(fontsize=10)),

                        # Heatmap width and height
                        # width = unit(16,"in"),
                        # height = unit(18,"in"),

                        # Column label parameters
                        column_labels = col_labels, column_names_rot = 45,

                        # Split columns by phase and row by activity type
                        column_split = factor(c(rep("Light",4), rep("Transition",3),rep("Dark",4),rep("Light+Dark",2)),
                                              levels=c("Light","Transition","Dark","Light+Dark")),
                        # row_split = factor(c(rep("1",5), rep("2",5))),

                        # Specify some column parameters.
                        column_title = NULL,
                        top_annotation = column_ha,

                        # Add signal direction arrows.
                        cell_fun = cell_fun1,

                        # Specify some parameters row dendrogram aesthetics and row labels.
                        row_title = NULL,
                        row_title_rot = 0,
                        #row_split = 3,
                        row_names_side = "left",

                        # Clustering parameters.
                        cluster_columns = FALSE,
                        cluster_rows = FALSE,

                        # Font sizes
                        row_names_gp = grid::gpar(fontsize=8),
                        column_names_gp = grid::gpar(fontsize=10),
                        row_title_gp = grid::gpar(fontsize=10),
                        column_title_gp = grid::gpar(fontsize=10)
                        )

ht <- draw(htlist2, merge_legend = FALSE, annotation_legend_list = lgd_list,
           annotation_legend_side = "top", align_annotation_legend = "heatmap_center")
ro <- row_order(ht)
co <- column_order(ht)

# Highlight startle endpoints and average acceleration in Dark endpoint
decorate_heatmap_body(paste0("BMC log(","\U03BC","M)"), row_slice = 1, column_slice = 1, {
  grid.rect(unit(1.025,'npc'), unit(0.7,'npc'),
            width = (length(co[[1]]) + length(co[[2]]))/length(co[[1]]) * unit(0.415,'npc') + unit(1,'mm'),
            height = (length(ro[[1]]) + length(ro[[2]]))/length(ro[[1]]) * unit(0.344,'npc') + unit(1,'mm'),
            gp = gpar(lwd=2, lty=1, fill=NA, col="red"), just=c('left','top')
  )
}) 
decorate_heatmap_body(paste0("BMC log(","\U03BC","M)"), row_slice = 1, column_slice = 1, {
  grid.rect(unit(2.05,'npc'), unit(1,'npc'),
            width = (length(co[[1]]) + length(co[[2]]))/length(co[[1]]) * unit(0.13,'npc') + unit(1,'mm'),
            height = (length(ro[[1]]) + length(ro[[2]]))/length(ro[[1]]) * unit(0.245,'npc') + unit(1,'mm'),
            gp = gpar(lwd=2, lty=1, fill=NA, col="red"), just=c('left','top')
  )
})

Supplemental Table 3: Number of hits per chemical.

count.hits.chm <- tcpl_out.dt[hitcall>0.8 & !endp%in%exclude, .N, by = .(name)]
datatable(count.hits.chm[order(-N)], colnames=c("Chemical Name","Number of Hits"), caption = "Number of Hits per Chemical",
          rownames = FALSE)

Supplemental Figure 5: Create plots for comparing Paraquat and Heptachlor epoxide activity profiles.

chemical <- "Heptachlor epoxide"

# Units
unit.t = "min"
unit.mov = "cm"
unit.conc = paste0("\U03BC","M")
prsp = "SA"
no.A = 10

# Extract chemical data
group <- unique(lmr0.egid[cpid == chemical, egid])

## Identify movement columns of interest
t.cols <- grep("vt", names(lmr0.egid), value = TRUE)
cols <- t.cols[(no.A+1):length(t.cols)]
A.cols <- t.cols[!(t.cols%in%cols)]

## extract data to be plotted, exclude acclimation
to.fit <- lmr0.egid[cpid==chemical | (wllt=="v" & egid==group), -A.cols, with=FALSE]

# create appropriate axes titles for plots
label.y <- "Speed"

# Format data for plotting

## calculate mean and 50% CIs for each vector column by concentration group, excluding concentration
exclude.A <- t.cols[!(t.cols%in%A.cols)]
means <- to.fit[, lapply(.SD, function(col) mean(col,na.rm=T)),
                .SDcols = exclude.A,
                by = conc]

## calculate CI's for transformed values then transform back
shift <- 1
logCIs <- to.fit[, lapply(.SD, function(x) log10(x+shift)), .SDcols=exclude.A, by=conc][
  , lapply(.SD, function(x) t.test(x,conf.level=0.50)$conf.int), .SDcols=exclude.A, by=conc]
CIs <- logCIs[, lapply(.SD, function(x) (10^x)-shift), by=conc][,lapply(.SD, function(col) abs(diff(col))/2), .SDcols=exclude.A, by=conc]

## elongate means and CIs data, and join
means_long <- data.table::melt(means, id.vars = "conc", variable.name = "t", value.name = "mean")
means_long[, t := sub("vt","",t)]
CIs_long <- data.table::melt(CIs, id.vars = "conc", variable.name = "t", value.name = "CI")
CIs_long[, t := sub("vt","",t)]
stats <- means_long[CIs_long, on = c("conc","t")][, conc := as.factor(conc)]
stats[, t := as.numeric(t)]

# create standard error of mean estimates by time period and plot as ribbons or error bars

# plot time-series data

## create title, x- and y-axis titles, legend labels, and legend title
title <- paste0("Sample Averaged Time-Series for ", chemical)
title.t <- paste0("Time (",unit.t,")")
title.mean <- paste0("Mean ", label.y, " (",unit.mov,"/",unit.t,")")
conc.n <- to.fit[wllq==1, .N, by=.(conc)][order(conc)]
legend.labels <- paste0(conc.n$conc, ", n=", conc.n$N)

title.legend <- paste0("Concentration (", unit.conc, ")")

## get better colors for plotting
N <- length(unique(to.fit[,conc]))
colors <- viridis::viridis(N)

## create x-axis breaks and labels
m <- as.integer(max(means_long[,t]))
x.breaks <- seq(from=no.A,to=m,by=10)
x.labels1 <- 2*seq(from=no.A,to=m,by=10)
x.labels2 <- x.labels1 - 2
x.labels <- paste(x.labels2, x.labels1, sep="-")

## plot
plot.HepEpox <- ggplot() +
          geom_point(data = stats, aes(x=t, y=mean, color=as.factor(conc))) +
          geom_line(data = stats, aes(x=t, y=mean, color=conc, group=conc)) +
          scale_x_continuous(breaks = x.breaks, labels = x.labels) +
          scale_color_manual(values = colors, labels=legend.labels) +
          geom_ribbon(data = stats,
                      aes(x=t, ymax=mean+CI, ymin=mean-CI, group=conc, fill=conc),
                      alpha = 0.25) +
          geom_rect(aes(xmin=10,xmax=30,ymin=-1.25,ymax=-0.5) ,fill="white", color="black") +
          geom_rect(aes(xmin=30,xmax=50,ymin=-1.25,ymax=-0.5) ,fill="black", color="black") +
          annotate("text", x=c(20,40), y=rep(-0.875,2), color=c("black","white"), label=c("Light Phase","Dark Phase"), size=3) +
          scale_fill_manual(values = colors, labels=legend.labels) +
          labs(title = title, subtitle = "Acclimation Period Excluded: 50% Confidence Bands",
               x = title.t, y = title.mean, color = title.legend) +
          guides(fill = "none") +
          theme_bw() +
          theme(text = element_text(size = 14))
signal.HepEpox <- tcpl_out.dt[name=="Heptachlor epoxide" & hitcall>0.8 & !endp%in%exclude, .(endp,top_over_cutoff)]
table.HepEpox <- gridExtra::tableGrob(signal.HepEpox, rows = NULL, cols = c("Endpoint Abbreviation", "Top Over Cutoff"))
chemical <- "Paraquat"

# Units
unit.t = "min"
unit.mov = "cm"
unit.conc = paste0("\U03BC","M")
prsp = "SA"
no.A = 10

# Extract chemical data
group <- unique(lmr0.egid[cpid == chemical, egid])

## Identify movement columns of interest
t.cols <- grep("vt", names(lmr0.egid), value = TRUE)
cols <- t.cols[(no.A+1):length(t.cols)]
A.cols <- t.cols[!(t.cols%in%cols)]

## extract data to be plotted, exclude acclimation
to.fit <- lmr0.egid[cpid==chemical | (wllt=="v" & egid==group), -A.cols, with=FALSE]

# create appropriate axes titles for plots
label.y <- "Speed"

# Format data for plotting

## calculate mean and 50% CIs for each vector column by concentration group, excluding concentration
exclude.A <- t.cols[!(t.cols%in%A.cols)]
means <- to.fit[, lapply(.SD, function(col) mean(col,na.rm=T)),
                .SDcols = exclude.A,
                by = conc]

## calculate CI's for transformed values then transform back
shift <- 1
logCIs <- to.fit[, lapply(.SD, function(x) log10(x+shift)), .SDcols=exclude.A, by=conc][
  , lapply(.SD, function(x) t.test(x,conf.level=0.50)$conf.int), .SDcols=exclude.A, by=conc]
CIs <- logCIs[, lapply(.SD, function(x) (10^x)-shift), by=conc][,lapply(.SD, function(col) abs(diff(col))/2), .SDcols=exclude.A, by=conc]

## elongate means and CIs data, and join
means_long <- data.table::melt(means, id.vars = "conc", variable.name = "t", value.name = "mean")
means_long[, t := sub("vt","",t)]
CIs_long <- data.table::melt(CIs, id.vars = "conc", variable.name = "t", value.name = "CI")
CIs_long[, t := sub("vt","",t)]
stats <- means_long[CIs_long, on = c("conc","t")][, conc := as.factor(conc)]
stats[, t := as.numeric(t)]

# create standard error of mean estimates by time period and plot as ribbons or error bars

# plot time-series data

## create title, x- and y-axis titles, legend labels, and legend title
title <- paste0("Sample Averaged Time-Series for ", chemical)
title.t <- paste0("Time (",unit.t,")")
title.mean <- paste0("Mean ", label.y, " (",unit.mov,"/",unit.t,")")
conc.n <- to.fit[wllq==1, .N, by=.(conc)][order(conc)]
legend.labels <- paste0(conc.n$conc, ", n=", conc.n$N)

title.legend <- paste0("Concentration (", unit.conc, ")")

## get better colors for plotting
N <- length(unique(to.fit[,conc]))
colors <- viridis::viridis(N)

## create x-axis breaks and labels
m <- as.integer(max(means_long[,t]))
x.breaks <- seq(from=no.A,to=m,by=10)
x.labels1 <- 2*seq(from=no.A,to=m,by=10)
x.labels2 <- x.labels1 - 2
x.labels <- paste(x.labels2, x.labels1, sep="-")

## plot
plot.Paraquat <- ggplot() +
          geom_point(data = stats, aes(x=t, y=mean, color=as.factor(conc))) +
          geom_line(data = stats, aes(x=t, y=mean, color=conc, group=conc)) +
          scale_x_continuous(breaks = x.breaks, labels = x.labels) +
          scale_color_manual(values = colors, labels=legend.labels) +
          geom_ribbon(data = stats,
                      aes(x=t, ymax=mean+CI, ymin=mean-CI, group=conc, fill=conc),
                      alpha = 0.25) +
          geom_rect(aes(xmin=10,xmax=30,ymin=-1.25,ymax=-0.5) ,fill="white", color="black") +
          geom_rect(aes(xmin=30,xmax=50,ymin=-1.25,ymax=-0.5) ,fill="black", color="black") +
          annotate("text", x=c(20,40), y=rep(-0.875,2), color=c("black","white"), label=c("Light Phase","Dark Phase"), size=3) +
          scale_fill_manual(values = colors, labels=legend.labels) +
          labs(title = title, subtitle = "Acclimation Period Excluded: 50% Confidence Bands",
               x = title.t, y = title.mean, color = title.legend) +
          guides(fill = "none") +
          theme_bw() +
          theme(text = element_text(size = 14))
signal.Paraquat <- tcpl_out.dt[name=="Paraquat" & hitcall>0.8 & !endp%in%exclude, .(endp,top_over_cutoff)]
table.Paraquat <- gridExtra::tableGrob(signal.Paraquat, rows = NULL, cols = c("Endpoint Abbreviation", "Top Over Cutoff"))
cowplot::plot_grid(plot.HepEpox, plot.Paraquat, table.HepEpox, table.Paraquat, ncol = 2, nrow = 2,rel_heights = c(3,1))

Potential Reference Chemicals

Supplemental Table 4: Active Chemicals Inducing the Strongest Signal in Endpoints

aenms <- strgEffectors_byEndp[, unique(aenm)]

# Format data and print as a good table
strgEffectors_formatted <- strgEffectors_byEndp[order(aenm), .(aenm,direction,cpid,conc,rval)]
strgEffectors_formatted[, rval:=round(rval, digits=2)]
setnames(strgEffectors_formatted, names(strgEffectors_formatted),
         c("Assay Endpoint Name","Signal Direction","Chemical Name",paste0("Concentration ","(\U03BC","M)"),"SEs from Control Mean"))
datatable(strgEffectors_formatted)

BMC Estimates

Table 4: Which chemicals appear most potent by median BMC across active endpoints?

select <- apply(bmd, 1, function(row) any(row<4))
med_BMC <- bmd[select,]
med_BMC[med_BMC==4] <- NA

med_BMC.1 <- apply(med_BMC, 1, median, na.rm=TRUE)
med_BMC.2 <- 10^(med_BMC.1[order(med_BMC.1)])

medBMC.dt <- data.table("Chemical Name" = names(med_BMC.2), "Median BMC" = signif(med_BMC.2, digits=3))

datatable(medBMC.dt, rownames = FALSE, 
          caption = "Median BMCs of Chemicals Across Active Endpoints",
          colnames = c("Chemical Name", paste0("Median BMC (","\U03BC","M)")))

Table 5: BMC ranges for chemicals.

BMDRanges[, BMCrange := signif(BMDRange, digits=3)]
datatable(BMDRanges, rownames = FALSE,
          caption = "BMC Ranges for Chemicals with Multiple Hits",
          colnames = c("Chemical Name", paste0("Range of BMCs (","\U03BC","M)")))

Supplemental Figure 6: Diazepam behavior data and active curve-fits.

plot_SA <- function(chemical) {
              # Units
              unit.t = "min"
              unit.mov = "cm"
              unit.conc = paste0("\U03BC","M")
              prsp = "SA"
              no.A = 10

              # Extract chemical data
              group <- unique(lmr0.egid[cpid == chemical, egid])

              ## Identify movement columns of interest
              t.cols <- grep("vt", names(lmr0.egid), value = TRUE)
              cols <- t.cols[(no.A+1):length(t.cols)]
              A.cols <- t.cols[!(t.cols%in%cols)]

              ## extract data to be plotted, exclude acclimation
              to.fit <- lmr0.egid[cpid==chemical | (wllt=="v" & egid==group), -A.cols, with=FALSE]

              # create appropriate axes titles for plots
              label.y <- "Speed"

              # Format data for plotting

              ## calculate mean and 50% CIs for each vector column by concentration group, excluding concentration
              exclude.A <- t.cols[!(t.cols%in%A.cols)]
              means <- to.fit[, lapply(.SD, function(col) mean(col,na.rm=T)),
                              .SDcols = exclude.A,
                              by = conc]

              ## calculate CI's for transformed values then transform back
              shift <- 1
              logCIs <- to.fit[, lapply(.SD, function(x) log10(x+shift)), .SDcols=exclude.A, by=conc][
                , lapply(.SD, function(x) t.test(x,conf.level=0.50)$conf.int), .SDcols=exclude.A, by=conc]
              CIs <- logCIs[, lapply(.SD, function(x) (10^x)-shift), by=conc][,lapply(.SD, function(col) abs(diff(col))/2), .SDcols=exclude.A, by=conc]

              ## elongate means and CIs data, and join
              means_long <- data.table::melt(means, id.vars = "conc", variable.name = "t", value.name = "mean")
              means_long[, t := sub("vt","",t)]
              CIs_long <- data.table::melt(CIs, id.vars = "conc", variable.name = "t", value.name = "CI")
              CIs_long[, t := sub("vt","",t)]
              stats <- means_long[CIs_long, on = c("conc","t")][, conc := as.factor(conc)]
              stats[, t := as.numeric(t)]

              # create standard error of mean estimates by time period and plot as ribbons or error bars

              # plot time-series data

              ## create title, x- and y-axis titles, legend labels, and legend title
              title <- paste0("Sample Averaged Time-Series for ", chemical)
              title.t <- paste0("Time (",unit.t,")")
              title.mean <- paste0("Mean ", label.y, " (",unit.mov,"/",unit.t,")")
              conc.n <- to.fit[wllq==1, .N, by=.(conc)][order(conc)]
              legend.labels <- paste0(conc.n$conc, ", n=", conc.n$N)

              title.legend <- paste0("Concentration (", unit.conc, ")")

              ## get better colors for plotting
              N <- length(unique(to.fit[,conc]))
              colors <- viridis::viridis(N)

              ## create x-axis breaks and labels
              m <- as.integer(max(means_long[,t]))
              x.breaks <- seq(from=no.A,to=m,by=10)
              x.labels1 <- 2*seq(from=no.A,to=m,by=10)
              x.labels2 <- x.labels1 - 2
              x.labels <- paste(x.labels2, x.labels1, sep="-")

              ## plot
              plot <- ggplot() +
                geom_point(data = stats, aes(x=t, y=mean, color=as.factor(conc))) +
                geom_line(data = stats, aes(x=t, y=mean, color=conc, group=conc)) +
                scale_x_continuous(breaks = x.breaks, labels = x.labels) +
                scale_color_manual(values = colors, labels=legend.labels) +
                geom_ribbon(data = stats,
                            aes(x=t, ymax=mean+CI, ymin=mean-CI, group=conc, fill=conc),
                            alpha = 0.25) +
                geom_rect(aes(xmin=10,xmax=30,ymin=-1.25,ymax=-0.5) ,fill="white", color="black") +
                geom_rect(aes(xmin=30,xmax=50,ymin=-1.25,ymax=-0.5) ,fill="black", color="black") +
                annotate("text", x=c(20,40), y=rep(-0.875,2), color=c("black","white"), label=c("Light Phase","Dark Phase"), size=5) +
                scale_fill_manual(values = colors, labels=legend.labels) +
                labs(title = title, subtitle = "Acclimation Period Excluded: 50% Confidence Bands",
                     x = title.t, y = title.mean, color = title.legend) +
                guides(fill = "none") +
                theme_bw() +
                theme(text = element_text(size = 18))

              return(plot)
}
chemical <- "Diazepam"
SAplot_Dzp <- plot_SA(chemical)
row1_Dzp <- rows_n[[chemical]][["avgS_L"]]
row2_Dzp <- rows_n[[chemical]][["RoA_L"]]
row3_Dzp <- rows_n[[chemical]][["avgS_T"]]

fit1_Dzp <- concRespCoreZR(row1_Dzp, do.plot = TRUE, verbose.plot = FALSE)[["plot"]]
fit2_Dzp <- concRespCoreZR(row2_Dzp, do.plot = TRUE, verbose.plot = FALSE)[["plot"]]
fit3_Dzp <- concRespCoreZR(row3_Dzp, do.plot = TRUE, verbose.plot = FALSE)[["plot"]]

p1 <- cowplot::plot_grid(fit1_Dzp, fit2_Dzp, fit3_Dzp, nrow=2)

cowplot::plot_grid(SAplot_Dzp, p1, nrow=2)

Figure 6: BMC confidence interval overlap.

tcpl_out.dt <- as.data.table(DNT60_tcpl_out)

to.plot <- tcpl_out.dt[hitcall>0.8 & endp%in%c("avgS_L","avgA_L","avgJ_L","hbt_L",
                                                 "strtlA","strtlAavg","strtlF",
                                                 "avgS_D","avgA_D","avgJ_D","hbt_D",
                                                 "avgS_T","AUC_r")]
to.plot[, endp := factor(endp, levels=rev(c("avgS_L","avgA_L","avgJ_L","hbt_L",
                                        "strtlA","strtlAavg","strtlF",
                                        "avgS_D","avgA_D","avgJ_D","hbt_D",
                                        "avgS_T","AUC_r")))]
to.plot[name=="Polybrominated diphenyl ether (PBDE)-47", name:="PBDE-47"]
ggplot(to.plot, aes(x=log10(bmd),y=endp)) +
  geom_point() +
  geom_errorbar(aes(xmin=log10(bmdl), xmax=log10(bmdu))) +
  scale_x_continuous(breaks=c(-3,-2,-1,0,1,2)) +
  facet_wrap(. ~ name)+
  theme_bw() +
  theme(text = element_text(size=14)) +
  labs(title = "Overlap of BMC Estimates for Chemicals", x = "log(BMC)", y = "Endpoint")

Supplmental Figure 7: D-Sorbitol active curve-fits.

chemical <- "Amphetamine"

row1 <- rows_n[[chemical]][["avgS_L"]]
row2 <- rows_n[[chemical]][["hbt1_L"]]
row3 <- rows_n[[chemical]][["strtlF"]]

fit1 <- concRespCoreZR(row1, do.plot = TRUE, verbose.plot = FALSE)[["plot"]]
fit2 <- concRespCoreZR(row2, do.plot = TRUE, verbose.plot = FALSE)[["plot"]]
fit3 <- concRespCoreZR(row3, do.plot = TRUE, verbose.plot = FALSE)[["plot"]]

legend <- cowplot::get_legend(fit1)
title <- cowplot::ggdraw() +
  cowplot::draw_label("Amphetamine Exposure: Active Curve-Fits",
                      x = 0, hjust = 0, size=24) +
  theme(plot.margin = margin(0,0,0,7))

edit_plot <- function(plot) {
  title <- plot$labels$title
  plot$labels$title <- gsub(paste0(chemical," for "), "", title)
  plot$labels$subtitle <- NULL
  plot <- plot + theme(legend.position="none")
  return(plot)
}

plots <- lapply(list(fit1, fit2, fit3), edit_plot)

p1 <- cowplot::plot_grid(plotlist=plots, nrow=2)
p2 <- cowplot::plot_grid(p1, legend, rel_widths = c(3, .4))

cowplot::plot_grid(title, p2, ncol=1, rel_heights = c(0.1,1))

Figure 7: BMC Heatmap

draw(htlist, merge_legend = FALSE, annotation_legend_list = lgd_list, annotation_legend_side = "top", align_annotation_legend = "heatmap_center")

Concordance Between this Study and Other Studies

How many chemicals had papers associated with them?

litReviewSummary[Number.of.associated.publications!=0, .N]
## [1] 37

How many papers of the 36 papers report chemical effect.

length( litReviewConc[neurotox==1 & DOI!="DNT60", unique(DOI)] )
## [1] 30

Supplmental Table 5: Number of Associated Chemicals and Reproducibility of Active Chemicals

litReviewSummary1 <- litReviewSummary[,-c(2,3,4)]
DT::datatable(litReviewSummary1, caption="Formatted Literature Review table.",
              extensions = 'FixedColumns',
              colnames = c("Chemical Name", "Number of Reviewed Publications", "Number of Reviewed Publications Reporting Effect", "Percent of Publications Reporting Effect"),
              options = list(scrollX=TRUE),
              rownames = FALSE)

Find median LOELs for DNT60 study.

DNT60_medLOELs <- DNT60LOELs[, .(cpid, med_LOEL = matrixStats::rowMedians(as.matrix(.SD),na.rm=TRUE)), .SDcols=names(DNT60LOELs)[-1]]

Table 6: Use median gabi LOEL to compare to other studies.

compareLOELs <-litReviewLOELs[DOI != "DNT60"][DNT60_medLOELs, on = .(cpid)]
setnames(compareLOELs, "med_LOEL", "DNT60_med_LOEL")
compareLOELs[, LOEL_diff := neuro_LOEL - DNT60_med_LOEL]

compareLOEL_summary <- compareLOELs[, .(mean_LOEL_diff = mean(LOEL_diff,na.rm=T), sd_diff_mag = sd(LOEL_diff,na.rm=TRUE)), by = .(cpid)][order(sd_diff_mag)]

compareLOEL_summary[, `:=` (mean_LOEL_diff = signif(mean_LOEL_diff,digits=3), sd_diff_mag = signif(sd_diff_mag,digits=3))]

DT::datatable(compareLOEL_summary,
              caption="Summarize Differences in LOEL Magnitudes",
              extensions = 'FixedColumns',
              options = list(scrollX=TRUE),
              rownames = FALSE)

The Effect of Developmental Toxicity on the Detection of Neurotoxicity

How many chemicals have at least one associated paper that reports chemical activity?

litReview_actives <- litReviewConc[DOI!="DNT60" & neurotox==1, unique(cpid)]
length(litReview_actives)
## [1] 31
# Gather data for later plots and questions
DNT60LOELs.melt <- melt(DNT60LOELs, variable.name = "endpoint", value.name = "LOEL")
DNT60LOELs.melt[, `:=` (DOI="DNT60", category="intra_neuro")]

compareLOELs <- litReviewLOELs[DOI != "DNT60",-"dev_LOEL", with=FALSE][, .(cpid, DOI, LOEL=neuro_LOEL, category="inter_neuro")]
dev_LOELs <- litReviewLOELs[, .(LOEL=unique(dev_LOEL),category="intra_dev"), by=.(cpid)]
compareLOELs <- rbind(DNT60LOELs.melt, compareLOELs, dev_LOELs, fill = TRUE)
gabi_LOEL_actives <- compareLOELs[category=="intra_neuro" & !is.na(LOEL), unique(cpid)]

discordant_inactives <- litReviewSummary[!cpid%in%gabi_LOEL_actives, .(cpid, Number.of.associated.publications, Number.of.publications.with.activity, `%.of.papers.with.activity`)][`%.of.papers.with.activity`>50][order(-`%.of.papers.with.activity`)]

How many chemicals with a gabi behavioral LOEL had an associated publication report observed effect at at least one concentration.

table(gabi_LOEL_actives %in% litReview_actives)
## 
## FALSE  TRUE 
##     6    15

What chemicals were active in other publications (>50%) but inactive in our analysis?

discordant_inactives
##                       cpid Number.of.associated.publications
##  1: Bis(tributyltin) oxide                                 2
##  2:       Cadmium chloride                                 1
##  3:               Caffeine                                 1
##  4:               Dieldrin                                 2
##  5:              Isoniazid                                 1
##  6:                  Maneb                                 1
##  7:               Nicotine                                 3
##  8:                 Phenol                                 1
##  9:           Lead acetate                                 5
## 10:             Acrylamide                                 3
## 11:             Colchicine                                 3
##     Number.of.publications.with.activity %.of.papers.with.activity
##  1:                                    2                 100.00000
##  2:                                    1                 100.00000
##  3:                                    1                 100.00000
##  4:                                    2                 100.00000
##  5:                                    1                 100.00000
##  6:                                    1                 100.00000
##  7:                                    3                 100.00000
##  8:                                    1                 100.00000
##  9:                                    4                  80.00000
## 10:                                    2                  66.66667
## 11:                                    2                  66.66667
discordant_inactives[Number.of.associated.publications > 1]
##                      cpid Number.of.associated.publications
## 1: Bis(tributyltin) oxide                                 2
## 2:               Dieldrin                                 2
## 3:               Nicotine                                 3
## 4:           Lead acetate                                 5
## 5:             Acrylamide                                 3
## 6:             Colchicine                                 3
##    Number.of.publications.with.activity %.of.papers.with.activity
## 1:                                    2                 100.00000
## 2:                                    2                 100.00000
## 3:                                    3                 100.00000
## 4:                                    4                  80.00000
## 5:                                    2                  66.66667
## 6:                                    2                  66.66667

Table 7: Median differences between developmental LOELs and study behavioral LOELs.

LOEL_diff_DN <- dev_LOELs[,.(cpid,dev_LOEL=LOEL)][compareLOELs[category == "inter_neuro", .(cpid, DOI, LOEL)], on = .(cpid)][, .(cpid, DOI, mag_diff=LOEL-dev_LOEL)]
devBhvLOEL.diff <- LOEL_diff_DN[, .(median_mag_diff = median(mag_diff,na.rm=T), sd_mag_diff = sd(mag_diff,na.rm=T)), by = .(cpid)][order(-median_mag_diff)]

DT::datatable(devBhvLOEL.diff[!is.na(median_mag_diff), .(cpid = cpid, median_mag_diff = signif(median_mag_diff,digits=3), sd_mag_diff = signif(sd_mag_diff,digits=3))], caption="Median Difference in Magnitude from Developmental LOEL to Literature Review Behavioral LOELs",
              colnames = c("Chemical Name", "Median Magnitude Difference", "SD in Magnitude Difference"),
              extensions = 'FixedColumns',
              options = list(scrollX=TRUE),
              rownames = FALSE)

Exclude those without majority of publications reporting activity.

cpid_maj_active <- litReviewSummary[`%.of.papers.with.activity`>50 & Number.of.associated.publications>1, cpid]
devBhvLOEL.diff[cpid %in% cpid_maj_active]
##                       cpid median_mag_diff sd_mag_diff
##  1:               Dieldrin       0.6365006   0.3373757
##  2:               Nicotine       0.3124595   0.1118929
##  3:             Heptachlor      -0.1020600   0.7071068
##  4: Bis(tributyltin) oxide      -0.2665894   1.0768379
##  5:           Lead acetate      -0.4586073   1.4744919
##  6:             Colchicine      -0.4900526   1.0098438
##  7:              Valproate      -0.6629268   0.3797085
##  8:            Haloperidol      -1.0401342   0.1938097
##  9:      Bisphenol A (BPA)      -1.4259687   0.5599005
## 10:               Diazepam      -1.6020600   0.0000000
## 11:           Chlorpyrifos      -1.6243364   1.1082443
## 12:             Fluoxetine      -1.9705933   0.5671248
## 13:          Acetaminophen              NA          NA
## 14:             Acrylamide              NA          NA

Figure 8: Developmental and nuero- toxic LOELs observed in this study and others.

compareLOELs1 <- data.table::copy(compareLOELs)

cpid_wO_pub <- litReviewSummary[Number.of.associated.publications==0, unique(cpid)]
cpid_w_pub <- litReviewSummary[Number.of.associated.publications!=0, unique(cpid)]

# Compare LOELs has different spellings
# Bis(tributyltin) oxide
# Diethylene glycol

compareLOELs1[cpid%in%cpid_w_pub, facet := 1]
compareLOELs1[cpid%in%cpid_wO_pub, facet := 2]

# Order chemicals in each facet first by developmental toxicity and then by median LOEL were devTox doesn't exist and then by median study LOEL where median LOEL doesn't exist
factor.orderA1 <- rev( compareLOELs[cpid%in%cpid_w_pub & category=="intra_dev" & !is.na(LOEL)][order(-LOEL), unique(cpid)] )
factor.orderB1 <- rev( compareLOELs[cpid%in%cpid_wO_pub & category=="intra_dev" & !is.na(LOEL)][order(-LOEL), unique(cpid)] )

factor.orderA2 <- DNT60_medLOELs[!(cpid%in%factor.orderA1) & cpid%in%cpid_w_pub & !is.na(med_LOEL)][order(med_LOEL), unique(cpid)]
factor.orderB2 <- DNT60_medLOELs[!(cpid%in%factor.orderB1) & cpid%in%cpid_wO_pub & !is.na(med_LOEL)][order(med_LOEL), unique(cpid)]

factor.orderA3 <- compareLOELs[!(cpid%in%factor.orderA2) & !(cpid%in%factor.orderA1) & cpid%in%cpid_w_pub & category=="inter_neuro", .(med_LOEL = median(LOEL,na.rm=TRUE)), by = .(cpid)][!is.na(med_LOEL)][order(med_LOEL), unique(cpid)]
factor.orderB3 <- compareLOELs[!(cpid%in%factor.orderB2) & !(cpid%in%factor.orderB1) & cpid%in%cpid_wO_pub][order(-cpid), unique(cpid)]

factor.orderA4 <- compareLOELs[!(cpid%in%factor.orderA3) & !(cpid%in%factor.orderA2) & !(cpid%in%factor.orderA1) & cpid%in%cpid_w_pub][order(-cpid), unique(cpid)]

factor.order <- c(factor.orderA4, factor.orderA3, factor.orderA2, factor.orderA1, factor.orderB3, factor.orderB2, factor.orderB1)

compareLOELs1$cpid <- factor(compareLOELs$cpid, levels=factor.order)
facet_labels = list("1" = "Chemicals with Associated Publications",
                   "2" = "Chemicals without Associated Publications")
facet_labeller <- function(variable, value) {
  return(facet_labels[value])
}
ggplot() +
  geom_point(compareLOELs1, mapping = aes(x=LOEL, y=cpid, color=category, shape=category), size = 3) +
  scale_shape_discrete(labels = c("Extra-Study Behavioral LOELs","Developmental LOEL","Study Behavioral LOELs")) +
  scale_color_manual(labels = c("Extra-Study Behavioral LOELs","Developmental LOEL","Study Behavioral LOELs"), values=colors[c(1,8,4)]) +
  theme_bw() +
  theme(text = element_text(size=14)) +
  facet_wrap(facet~., dir = "v", as.table = TRUE, drop = TRUE, scale = "free_y", strip.position = "right", labeller = facet_labeller) + 
  labs(shape = "Critical Doses",
       color = "Critical Doses",
       x = paste0("Concentration (log10(",paste0("\U03BC","M"),"))"),
       y = "Chemical Name")

Supplemtnal Figure 8: Unity plot of median potencies.

# plot medians with quantiles and interquartile range
cHairPlot.data1 <- compareLOELs[category=="intra_neuro", .(med=summary(LOEL)[3], Q1=summary(LOEL)[1], Q3=summary(LOEL)[5]), by=.(cpid)]
cHairPlot.data2 <- compareLOELs[category=="inter_neuro", .(med_x=summary(LOEL)[3], Q1_x=summary(LOEL)[1], Q3_x=summary(LOEL)[5]), by=.(cpid)]
# _x indicates values are derived from "extra" laboratory studies

cHairPlot.data <- merge(cHairPlot.data1, cHairPlot.data2, all = TRUE)
ggplot(cHairPlot.data, aes(x=med, y=med_x)) +
  geom_point() + 
  geom_abline(intercept=c(-1,0,1), slope = c(1,1,1), linetype=c("dashed","solid","dashed")) +
  ggrepel::geom_text_repel(aes(label=cpid,color=cpid), box.padding = 0.5) +
  lims(x=c(-2,2), y=c(-2,2)) +
  labs(x = "Median Neurotoxic LOEL Observed in this Study",
       y = "Median Neurotoxic LOEL Derived from Literature Review Studies",
       title = "Concordance of Neurotoxic LOELs Across Studies") +
  theme_bw() +
  theme(legend.position="none")